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BIOLOGICAL AGGREGATION DRIVEN BY SOCIAL AND 
ENVIRONMENTAL FACTORS: A NONLOCAL MODEL AND ITS 
DEGENERATE CAHN-HILLIARD APPROXIMATION 

ANDREW J. BERNOFF t§ AND CHAD M. TOPAZ* 

Abstract. Biological aggregations such as insect swarms and bird flocks may arise from a com¬ 
bination of social interactions and environmental cues. We focus on nonlocal continuum equations, 
which are often used to model aggregations, and yet which pose significant analytical and computa¬ 
tional challenges. Beginning with a particular nonlocal aggregation model [Topaz et ai, Bull. Math. 
Bio., 2006], we derive the minimal well-posed long-wave approximation, which is a degenerate Cahn- 
Hilliard equation. Energy minimizers of this reduced, local model retain many salient features of 
those of the nonlocal model, especially for large populations and away from an aggregation’s bound¬ 
aries. Using the Cahn-Hilliard model as a testbed, we investigate the degree to which an external 
potential modeling food sources can be used to suppress peak population density, which is essential 
for controlling locust outbreaks. A random distribution of food sources tends to increase peak density 
above its intrinsic value, while a periodic pattern of food sources can decrease it. 

Key words, partial integrodifferential equation, Cahn-Hilliard, biological aggregation, energy, 
minimizer, locust 


1. Introduction. This paper has three primary aims. First, beginning with a 
partial integrodifferential equation model of biological aggregation, we show that a 
long-wave approximation yields a degenerate Cahn-Hilliard equation akin to models 
describing phase separation in materials science and evolution of thin fluid films. 
Second, we study solutions of the degenerate Cahn-Hillard model and demonstrate 
that they compare well with those of the original model. Using the Cahn-Hilliard 
model is advantageous because it eliminates many of the analytical and computational 
challenges of nonlocal equations, particularly in higher dimensions. Our methodology 
is to study the energy minimizers of this local equation, rather than studying the 
time evolution of the original nonlocal equation. Finally, we use the Cahn-Hilliard 
model as a testbed to assess whether imposing an external potential modeling the 
environment, e.g., food sources, can be used to control peak density in aggregations, 
which is of interest for locust swarms. 

Biological aggregations such as bird flocks, fish schools, and insect swarms are 
driven by social interactions between group members. Typical social forces include 
attraction, repulsion, and alignment [3311371 which are preprogrammed by evolution 
and which are facilitated by sight, sound, smell, touch, or some combination of these 
senses. The accepted biological paradigm is that repulsion operates over shorter dis¬ 
tances than attraction, but that it is much stronger over those distances. In the 
context of many mathematical models, these are necessary conditions to achieve co¬ 
hesive groups with finite densities [511 [sg. 

Because sensing operates over finite rather than infinitesimal distances, it is most 
naturally described via nonlocal operators. The most common mathematical descrip¬ 
tions assume that social interactions take place in an additive pairwise manner, and 
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that the strength of interaction of organisms depends on the distance between them. 
We will use a kinematic model, where an individual’s velocity is determined by so¬ 
cial interactions dependent on the positions of other organisms; we will eventually 
augment this model to allow taxis due to environmental effects. 

For a population that is well-described by a continuum density u(x, t), the social 
velocity v, can be modeled as 

v(x,t) = y f(x- y)'u(y,t)dy = f * u. (1) 

Because the interactions are pairwise, additive and solely a function of position, the 
velocity can be written as the convolution of a social kernel, f, that encodes the 
influence of organisms at location y on those at location x. Such convolution-type 
continuum aggregation models date back to the seminal work of [49j , which examines 

ft -I- V • (wv) = V = f * u, (2) 

in one spatial dimension. The left hand side of this equation is the material derivative 
of the population density u, that is, the rate of change of the population density 
moving at the velocity arising from social interactions. On the right hand side, D 
is a diffusion constant. The diffusion term was intended to model social repulsion 
operating over short distances, driving flux down population gradients. This model 
conserves mass, preserves positivity of the density, and can display an instability from 
a homogeneous state to patterned states. However, a key feature of biological swarms 
is that they are compactly supported with sharp edges, and cannot reproduce this 
property. This is due to the linear diffusion; even compactly supported initial state 
are instantaneously smoothed and densities become positive throughout the domain. 

More recently, a population model that has received much attention is the aggre¬ 
gation equation. 


u -f V • (uv) =0, V = — VQ * M, (3) 

where the velocity, v, can incorporate both attractive and repulsive social forces. The 
social kernel is written as minus the gradient of a social potential, —VQ. One typically 
assumes that the social potential is directionally isotropic and symmetric (the social 
force between two organisms is equal in magnitude and opposite in orientation) which 
allows us to write the social potential as Q = Q(r), solely a function of the distance 
between organisms, r. Typically Q{r) is decreasing for small r to model repulsion, 
increasing over larger r to model attraction, and asymptotically flat for large r as 
organisms cannot sense each other at large distances. 

The aggregation equation has recently been the subject of a rich and rapidly 
growing literature, including [TUI HU [131 [TSHT51 HH |331 IMl SOI H5M51 1531 IBS] and 
has also been adapted to describe specific biological organisms, such as locusts [BB] . 
Crucially, ([^ minimizes an energy, 

Eiu) = ^J y M(x)Q(|x-y|)M(y)dydx, (4) 

which can be used to analyze the existence and stability of equilibria [HI El [SOl |39l |BI] . 
Despite its popularity in the literature, one shortfall of this model is that it does not 
produce so-called well-spaced groups. In well-spaced biological swarms, individuals 
tend to pack only up to a maximum density, and no further. In contrast, for (§, if 
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there is a minimizer for a particular mass M, doubling the mass yields a minimizer 
with double the density. This is unbiological, and results from the attraction and the 
repulsion scaling equally with density. 

An alternative aggregation model that does form well-spaced groups is 

It -I- V • (uv) =0, V = —VQa * u — mVm, (5) 

initially proposed in [S^. While ^ models attraction and repulsion nonlocally, ([^ 
models attraction nonlocally via a kernel Qa, but models short-range repulsion via the 
(local) nonlinear diffusive velocity —uVu. This model is similar to ([^ but incorporates 
diffusion that is density-dependent and degenerate. That is, the diffusive flux tends 
to zero as the density tends to zero, which allows for compactly supported solutions. 

Eq. ([ 5 ]), like (§, allows the formation of compactly supported solutions 11121 Eg. 
It minimizes the energy, 

E{u) = ^ j J u(x)Qa(|x-y|)u(y)dydx-k ^ J dx. (6) 

The nonlinear diffusion has two notable effects. First, because the repulsion domi¬ 
nates at higher densities, the density is bounded as the mass M increases, leading to 
aggregations of approximately constant density at large mass [2l |65] . Second, the 
introduction of the gradient as part of the repulsion in the social velocity smooths the 
density profile at the support’s boundary. Solutions are continuous with a discountin- 
uous gradient, often referred to as a non-zero contact angle. Further work on ([^, and 
on a generalization incorporating repulsion of the form —has examined exis¬ 
tence and uniqueness maEn, regularity [32], global attractors |4T|, and properties 
of steady states as a function of the power p of the nonlinear repulsion [24j . 

Despite the promise of (§ as a model, it poses significant challenges. Beyond 
the analytical difficulties associated with the nonlocal operator, even numerical sim¬ 
ulation is problematic. First, for 2-d simulation on an n x n grid, the convolution 
requires 0{n^) operations for quadrature or 0{n'^ logn) for pseudospectral evaluation. 
Second, the steady states of ([^ can have steep edges, like some steady states of ([^. 
These require a high degree of spatial resolution. Third, standard methods produce 
oscillations emanating from the contact points, and such artifacts must be eliminated. 
Fourth, and relatedly, a numerical scheme must preserve the positivity of the model. 
It is notable but perhaps not surprising that fully two dimensional simulations of ([^ 
are rare in the literature |29[ I65j . At present, most published simulations are either in 
1-d or in a radially symmetric 2-d geometry. This is true of as well, excepting spe¬ 
cial cases where the computation reduces to a boundary integral problem describing 
the motion of a self-deforming curve [niiMj- Given the attractive features of ([^ from 
a modeling perspective, and yet given the analytical and computational challenges, 
we investigate whether there exists a simpler, local model that nonetheless preserves 
important features of ©■ 

In summary of this discussion, © produces groups whose density is bounded 
as mass M increases but the model is difficult to analyze and simulate. There is 
a long mathematical history of approximating dynamics in the limit of slow varia¬ 
tion to simplify governing equations |35[ I59j . As we will describe in detail in this 
paper, an expansion of the convolution in ([^ the long-wave limit yields, after non- 
dimensionalization. 


Qa*U 


—u — V^u. 


(7) 


4 


A.J. BERNOFF AND C.M. TOPAZ 


Retaining only the first term in this expansion yields a nonlinear diffusion equation 

tt + V • (itv) = 0, V = (1 — m)Vu, (8 ) 

which is ill-posed due to negative (backward) diffusion at small densities. However, 
retaining the second term in the expansion yields stabilizing fourth-order degenerate 
diffusion. We will concentrate on this degenerate fourth-order equation, 

ft H-V • (uv) = 0, -V = —uVu + V{u + V'^u). (9) 


It is curious that for this model, well-posedness at short wavelengths comes from the 
fourth-order term derived from nonlocal attraction. As attraction operates at longer 
wavelengths than repulsion, one must retain higher derivatives to model it faithfully. 
Mathematically, the fourth-order term describes attractive forces that are responsive 
not only to variations in density but to the convexity of the density. This yields a 
social force analogous to a surface tension, and has the effect of damping variation in 
the density. A moment expansion similar to the approximation Q is used to derive 
a different model for the movement of ecological populations in [52]. 

A great deal is known about As we will discuss below, it is related to both the 
Cahn-Hilliard equation [56] which arises as a model of phase separation and to thin- 
film equations which model surface-tension driven wetting of a substrate by a viscous 
fluid 155115^ . While the existence theory for our model is more strongly related 
to previous work on thin film equations, the dynamics are akin to phase separation 
and therefore we refer to our model as a Cahn-Hilliard model. Our model § also 
may be considered a degenerate Sivashinsky equation (57] |62| . 

After nondimensionalization, Q is equivalent to 


■u = V • 



( 10 ) 


which reveals a little more of its structure. It is the composition of a negative-definite 
second-order self-adjoint operator (V • uV) with the first variation of an energy, 

E{u) = j -|V«p-- + -dx. (11) 

This equation is an example of the more general class of Cahn-Hilliard models [56] 
which take the form, 


u = V • (M(u)V [f{u) - V^m] ) . 


( 12 ) 


Here, M{u) is a non-negative mobility, so V-M{u)V is a negative-definite second-order 
self-adjoint operator. The equation minimizes an energy 


E[n) = y'l|Vup + A(u) 


dx, 


(13) 


where F{u) = f f{u) du is most commonly taken to be a double well potential in the 
Cahn-Hilliard literature. Eq. (12) was first developed in [25l |28] to describe binary 


alloys. A hallmark of this model is that it exhibits phase separation, where the density 
segregates into regions in which u assumes one of the two minimum values of F, and 
these regions are separated by narrow transition layers [IHl [17| |33| . 
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Our model (10) borrows a feature more familiar in the thin him context. Bio¬ 
logically, the density u must remain nonnegative, much like a huid him thickness. In 
our model, this constraint is enforced by the degenerate mobility M{u) = u. This 
feature not only ensures that the density remains nonnegative, but it also allows re¬ 
gions of zero density, akin to dry spots in huid hlms. Degenerate mobilities have been 
considered in both the Cahn-Hilliard context [38l |56] and in the thin hlms context 

[siiTiiiiiiiiiiniiiniiiiiisiisn]. 

The work in [58] nicely summarizes and extends the analytical progress that has 
been made over the past two decades on a class of thin him equations including 
(10), primarily in one spatial dimension. For our model in one dimension, [58] shows 
that there exist nonnegative solutions with compact and potentially disjoint support. 
Where u is positive, these solutions are inhnitely differentiable. Moreover, these 
solutions are continuously differentiable everywhere. Stated differently, the slope at 
the contact point, where u touches down to zero, is zero. Additionally, the touchdown 
is generically quadratic and the speed of propagation of the contact points is hnite. To 
our knowledge, theory for uniqueness and existence in higher dimensions is lacking. 
Still, motivated by these one dimensional results, we will primarily consider solutions 
that are inhnitely diherentiable for positive values and touch down with zero contact 
slope. 

The rest of this paper is organized as follows. In Section we perform a long¬ 
wave expansion of ^ leading to the degenerate Cahn-Hilliard equation ([^ , augment¬ 
ing both models with an external potential describing the environment. Section [^ 
calculates the linear stability of the two models, builds a framework for examining 
them as energy minimization problems, and describes how numerical computation of 
minimizers will be performed. Section [^ studies in detail the minimizers for the case 
of no external potential. In brief, minimizers of both models share many properties, 
including compact support and a common peak density in the large mass limit. With 
the Cahn-Hilliard model established as a reasonable approximation of the nonlocal 
one, in Section [^ we consider an external potential and show that when a potential 
well is sufficiently steep, it can be replaced with an approximate obstacle potential 
consisting of regions where the potential is zero separated by inhibited regions where 
the potential is effectively infinite. Finally, Section uses these obstacle potentials 
to investigate how peak population density is affected by resource distribution. A 
random distribution of resources tends to increase peak density above its intrinsic 
value while a periodic distribution can decrease it. 


2. Derivation of the local model. We begin with the nonlocal model of [55] 
augmented with an external potential W (x) which models taxis due to environmental 
factors. 


■u -I- V • (uv) =0, V = —VQa * u — RvS/u — VW, (14) 

where u{'x.,t) is the population density and x G K”. The parameter i? > 0 controls 
the strength of short-range repulsion. For now, we assume the potential W to be 
continuously differentiable, although eventually we will consider a discontinuous limit. 
The kernel Qa in the convolution above describes social attraction over hnite (rather 
than inhnitesimal) distances. We place several assumptions on this kernel: 

1. Qa is radial, that is, Qa = Qa(r), f = |x|. This is the case for social interac¬ 
tions that are isotropic in space. 

2 . Q'a{r) > 0 so that Qa describes (pure) attraction. 
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3. As |r| —)■ oo, Qa approaches a constant since organisms at very long distances 
do not interact with each other. Without loss of generality, we assume Qa —>■ 0 
as |r| —)■ oo. 

4. Qa < 0 as a result of the previous two assumptions. 

This is a dimensioned model; by rescaling x, t, and u, we may write a dimensionless 
version of (14) as 


■u + V • (uv) =0, V = —\7Qa * u — vS/u — VW, (15) 

where we have chosen i? = 1. Also, we choose our nondimensionalization such that 
Qa has unit volume, that is, 

Mo= f Qa{r)d:si = -1, (16) 


where Mq denotes the zeroth moment. In addition, we specify the second moment of 
the potential in order that Qa has a length scale of order unity. This condition is 


M 2 = [ QaiQr"^ dx =-2n. 


(17) 


This rescaling, which we take without loss of generality, will be convenient for our 
subsequent calculations. 

Presently, we will approximate (151 with a local partial differential equation by 
applying the convolution theorem to — VQa *u and performing a long-wave expansion 
in Fourier space EH [Ml- Define the Fourier transform of a function /(x). 


Sif) = f= /We 


—/kx 


dx. 


(18) 


Since our kernel Qa is radial in physical space, its Fourier transform is radial in the 
transformed variable k. So 


Qa{k)= [ Qa(r)e-**^"dx, 


(19) 


where A: = |k| and r = |x|. Expanding the complex exponential, 

OO 

= E(-') 


zkx 


m—0 


(k-x)"" , 1,, ,9 

-;— = 1 — ik • X — - (k • x)^ 

ml 2^ ' 


and substituting into (19) will yield a moment expansion for Qa- 


the first three terms of this expansion. For the first term, 


[ Qair)dx = Mq =-1, 

jR" 


due to our choice of scaling. For the second term. 


-if (k • x)Qa(?') dx = 0, 

dR" 


+ • • • , ( 20 ) 

We will retain only 

( 21 ) 

(22) 
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because the term linear in k vanishes due to symmetry; in fact, all terms odd in k 
vanish. Finally, for the third (quadratic in k) term, we have 


1 

2 


^ n n r « 

{k-x)'^Qair)dx=--'^'^kpkq J Qa{r)xpXqdx 




-Ill'll 

p^l 

n 

p^l 

h‘^ 

= k\ 


Qa{r)xldyi 


Qa{r)r'^ dx 


(23a) 

(23b) 

(23c) 

(23d) 

(23e) 


Here, the right hand side of (23a) expands the square of the dot product as a double 
sum, (23b I utilizes symmetry to recognize the vanishing non-diagonal moments, (23c I 
uses radial symmetry of Qa, (23d) substitutes the definition of M 2 , and (23e) uses 
the normalization of M 2 in (17). 

Putting these results together. 


Qa = -l + k^ + 0{k^), (24) 

and therefore, we approximate 

Qa*u~ —u — V^u, (25) 


which constitutes a long-wave approximation. With this approximation, the original 

(26) 


governing equation (15) becomes the new, local equation 

li -I- V • (uv) =0, V = V(u -I- — W) — uVu. 


This equation may be rearranged as 


M = V • uV 


+ W — u — 


(27) 


which is a Cahn-Hilliard equation of type with degenerate mobility M(u) = it, 
f{u) = V?/2 — u, and an external potential bF(x). We will explore the degree to 
which minimizers of this truncated model approximate those of the nonlocal model 
( [Ts] ). While the former are easier to compute owing to the purely local interaction 
rules, it is important to realize what is lost in the local approximation. In (15), 


disjoint clumps of density interact with each other exponentially weakly owing to the 


nonlocality. However, in (27), disjoint clumps do not interact at all owing to the local 


behavior and degenerate diffusion. 

3. Basic model characteristics. To recapitulate, we are studying the nonlocal 
model 


ft = V • 



Qa * U -\- 



(28) 
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and the degenerate Cahn-Hilliard equation which is a long-wave approximation to it, 

1 


M = V • mV 


—M — V^M -|- -M^ -I- W 


(29) 


where M(x,t) is a non-negative density. To complete the formulation of the problem, 
one must specify the domain, the boundary conditions, the external potential W and 
initial conditions for the density m. For much of this paper we consider all of M" or 
periodic domains in one and two dimensions. Solutions to problems with degenerate 
diffusion often have compact support O HU |M] • For the Cahn-Hilliard problem ( [^ 
we expect u{x,t) to be continuously differentiable everywhere and infinitely differen¬ 
tiable for all X in the support of u. For the nonlocal model (28), we expect u{x,t) 
may have jump discontinuities at the edge of its support, but will be infinitely dif¬ 
ferentiable for all X in the support of u [HIS]. Leveraging the compact support of 
the solution, when we wish to consider densities confined to a finite domain rather 
than specifying boundary conditions at the edge of a domain, we will either consider 
a large periodic domain or we will specify that W (x) is sufficiently large to drive the 
density to zero except in some subdomain of interest. 

Both models are conservation laws which are well known to conserve total mass 
(see, e.g., [65]). If 


M = u dx. 


(30) 


where is the domain (or the support of the solution) then 

dM 


dt 


= 0 . 


(31) 


For many problems with a reflection symmetry, including our models (281 and 


(29), the location of the center of mass will also be conserved when the domain is 


and the mass is finite. In this case, for both (28) and (29), if we define 


1 f 

X = — / XM dx, 

M J a 


a straightforward calculation shows that 


dX 1 /• If 

—— = — / XMdx = —— / uVW dx, 
dt M M 


(32) 


(33) 


from which we deduce that in the absence of an external potential, the center of mass 
will be stationary. In the presence of an external potential, mass tends to migrate 
towards minima of the potential. 

3.1. Linear stability. In the absence of an external potential {W = 0), both 


(28) and (29) admit steady states with constant density m > 0. To analyze their linear 
stability, let m(x, t) = m -|- u(x,t), linearize the equations in m, and let u oc 
{k = |k|) to yield the dispersion relation a{k). For both models, (t( 0) = 0 for all u 
because the models are of conservation form 1351 . 


For the nonlocal model (28), 


cr(fc) = —k^u u — Q, 


(34) 
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a result also stated in [65]. The Fourier transform of the attractive kernel Qa depends 
on spatial dimension. From our nondimensionalization in Section 

Q, = -l + fc2 + 0(fc4), (35) 


This guarantees that the dispersion relations for (28) and (29) are identical to 

To further analyze the model one must specify Qa- A commonly used example 
which we also adopt is the Laplace potential. In our nondimensionalization, this 
potential is 




(36) 


where 


in which case 


_ ^/^T /n + 1\ 

” fW) 




n + 1 
2 


Qa{k) 



rr+1 

2 


(37) 


(38) 


which is a power of the Lorentzian function. We see that Qa(0) = — 1 and Qa{k) is 
increasing in fc, and tends to zero. Consequently, if u > 1, cr(/c) < 0 for all fc > 0 and 
the homogeneous steady state is linearly stable. However, if u < 1, there is a long¬ 
wave instability. There exists a hnite, continuous band of wave numbers extending 
from A: = 0 to 


k = 



(y-2/(n-M) _ 


(39) 


for which cr(fc) > 0. 

For the Cahn-Hilliard model (29), the calculations are signihcantly simpler. We 
obtain 


cr(fc) = —k'^u [(u — 1) -f /c^] . 


(40) 


The stability theory is analogous to the nonlinear model. If u > 1, cr(fc) < 0 for all 
/c > 0 and the homogeneous steady state is linearly stable. However, if u < I, there 
is a long-wave instability. There exists a finite, continuous band of wave numbers 
extending from k = 0 to k = y/l — u for which a(k) > 0. 

It is worthwhile to consider why a fourth-order truncation of (28) is superior to 
a second-order one. For a second-order truncation, the dispersion relation is simply 


a{k) = —k^u{u — I). 


(41) 


The stability threshold is still u = 1 as above. However, in the linearly unstable 
regime, all modes have positive growth rates which grow unboundedly with k, sim¬ 
ilar to the backwards heat equation. Thus, the linear problem for the second-order 
truncation is ill-posed in this case, shedding light on why an approximation to fourth 
order is desirable. Our fourth-order truncation is analogous to the strategy used in 
deriving amplitude equations for marginally long-wave unstable systems via modula¬ 
tion theory [35] . Retaining the destabilizing second-order and stabilizing fourth-order 
terms yields the most parsimonious truncation that is linearly well-posed. 
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3.2. Energy. We will show that both the nonlocal model (28) and the Cahn- 
Hillard model ( p^ can be interpreted as gradient flows. Both may be written in the 
form 


= V • uV [£« + f{u) + W], 


(42) 


where £ is a linear operator. For ([^), Cu = Qa * u and for (29), Cu = —— u. 
For both models, f{u) = u^/2. In (42), V • uV is a negative-definite self-adjoint 
second-order operator on non-negative functions u. 

We now define quantities useful for constructing an energy. First, define a sym¬ 
metric quadratic form for which £u is the first variation. For (28) define 

Q[u,v] = l^ [ uQa*vdx=l^ [ / (5a(|x-y|)M(x)i;(y)dxdy, (43) 

^ Jn ^ Jn Jn 


and for (29) define 


Next, choose 


Q[u, u] = ^ f \7u ■ \7v — uv dx. 

2 Jq 

F{u) = ( u^/6dx, 

Jn 


(44) 


(45) 


whose first variation is f{u). 

The total energy for both models is 


E{u) = Q[u,u] + F{u) + / Wudx. 

Jn 


The first variation of the energy is 


SE 

— ^ Cu + f{u) + W, 


and both evolution equations can be written as 

Ut = '^ ■ uV 


SE 

5u 


(46) 


(47) 


(48) 


It follows that the evolution is a gradient flow in a weighted metric [iisiisniiiiiniiisz]- 
The time derivative of the energy is 


dE 

dt 


In 


-^ut ax, 
ou 


f SE^ ^SE , 

/ -^V • mV dx. 
In aM Su 

SE ^ 


u 


m 


Su 


dx. 


(49a) 

(49b) 

(49c) 


Consequently, energy is decreasing except for stationary states where SE/Su is con¬ 
stant on the support of the solution. 
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3.3. Energy minimization and nnmerics. Because our governing equations 
are gradient flows, we expect almost all initial conditions to evolve to energy minimiz- 
ers. We exploit energy minimization to explore the equilibria of both models, namely 


(28) with the Laplace potential (36), and (29). This strategy plays a prominent role 
for the remainder of this paper. In order for a solution u to be an equilibrium, it 


must minimize the energy (46), satisfy the appropriate mass constraint, and be non¬ 


negative. In summary, the constrained minimization problem is 

minE{u) = Q[u,u\ + F{u) + / Wudx, 
“(x) Jn 

subject to the constraints: 


(50a) 


1 dx = M, 


u(x) > 0 for all x G ft. 


(50b) 

(50c) 


We solve this problem using the fmincon minimization routine of Matlab’s opti¬ 
mization toolbox, which accepts an objective function as well as constraints. We use 
periodic domains in one and two dimensions with equally spaced grid points. For the 
nonlocal model on a periodic domain, we replace the attractive potential Qa with its 
periodized analog and evaluate the convolution in Q via matrix multiplication. For 
the Cahn-Hilliard model, we evaluate derivatives using a centered difference on each 
interval in one dimension, or each cell in two dimensions. To aid convergence of the 
algorithm, we also supply the minimization algorithm with the gradient of our col¬ 
located approximation of the energy. We typically choose a random initial condition 
to seed the algorithm. Sometimes we solve the minimization problem on a coarse 
grid and then refine the grid by factors of two in order to obtain highly-resolved final 
states. 

Our methodology neglects the dynamics of the approach to equilibria, and relat- 
edly, the issue of whether those equilibria are accessible on relevant biological time 
scales. Nonetheless, we see the exploration of minimizers as an important first step 


in understanding (28) and (29). 


4. Energy minimizers in the absence of an external potential. 

4.1. Motivation and numerical computations. Recalling that the nonlocal 


equation (28) and the degenerate Cahn-Hilliard equation (29) have steady, compactly 


supported solutions and motivated by the aggregation of biological organisms, we now 
study the simplest case of a single clump having total mass M on an effectively infinite 


domain with no external potential {W = 0). A phase plane analysis of (28) suggests 


that in one spatial dimension, for each mass M, there is a single energy-minimizing 
clump [5S] . For a generalization of (28) to diffusion of power-law form, rigorous proofs 


show existence and uniqueness in higher dimensions, subject to appropriate assump¬ 


tions [SIS]. Similarly, degenerate Cahn-Hilliard models such as ( [2^ are gradient flows 
that minimize an energy; see [501156] and references therein. Our primary strategy 


for studying (28) and (29) will be to exploit energy minimization. 


To begin, we numerically compute minimizers as described in Section varying 
total mass M as our parameter for both models, and in both one and two dimen¬ 
sions. When varying M, we choose values of domain length L such that the spatially 
homogeneous solution is linearly unstable to a single mode. Strictly speaking, the 
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Fig. 1. Numerically calculated energy minimizers in one dimension. The nonlocal model 
result for (28^ appears as a dashed curve, and the Cahn-Hilliard model result for 1(29}^ appears as a 
solid curve, (a) Mass M = A, periodic domain length L = 40, calculated with N = 400 grid points. 
The minimizers have a well-defined peak, (b) M = 35, L = 40, N = 400. The minimizers have a 
plateau of peak density ||ti||oo ~ 1.5, and the two models coincide well except in the transition region. 


minimizers of the nonlocal model (281 are affected by their periodic images in our 
computation. In practice, with our choice of attractive potential (361, these effects 


are exponentially small and invisible on the plots shown. For the Cahn-Hilliard model 
(29), the compactly supported solution does not interact with its periodic images. 


For both models in one dimension, example minimizers appear in Figure 
Panel (a) compares small mass solutions and (b) compares large mass solutions. The 
small M solutions have well defined peaks, whereas the large M solutions acquire a 
plateau-like structure; they are all compactly supported. The computational results 


for (28) and (29) differ at small masses, but coincide for large M, except in the tran¬ 


sition region. Figure scans two key properties of the clumps as a function of M, 
namely the peak density ||u||oo and the size of the support || supp(u)||. The peak den¬ 
sity increases but saturates to a value of 1.5 with increasing M; we will later analyze 
this important property. The support increases with M, and the growth appears to 
be linear at large M (as we later confirm). 

Example minimizers in a fully two dimensional geometry appear in Figure The 
minimizers are radially symmetric. Panels (a) and (b) show clumps for small mass 
M for and ( [2^ respectively. To give a better sense of the mass distribution, 
Figure |4|(a) shows the radial profiles for ( [^ (dashed) and (29) (solid). Figures [^c) 
and (d) and Figure|^b) are analogous, but for large M. The results are similar to the 
one dimensional case. At small M, the clumps have a well-defined peak and at large 
M, they have a plateau of density ||u||oo ~ 1-5. The two models agree well, especially 
for large M. 

To better understand these results, it is helpful to analyze the steady state prob¬ 
lem. To find conditions for energy minimizers, we recall the discussion of Eli- Let 


i(x) = u-\- Su. 


(51) 


Here, u is an equilibrium solution of fixed mass M and 5u is a small perturbation of 
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M 


Fig. 2. Properties of energy minimizers of the nonlocal model {2^ (shown as squares) and the 
Cahn-Hilliard model ^29^ (shown as circles). For each mass M, and for both models, we calculate 
a minimizer on a periodic domain of length L (increasing with M so as to be effectively infinite) 
using lOL grid points. The dashed and solid curves are to guide the eye. The two models agree well 
for large M. (a) Peak density ||w||oo increases with M and saturates to a constant value of 1.5. 
(b) The size of the support, || supp(?i)||, increases and grows linearly at large M. 


zero mass, so that overall mass is conserved. Thus, 


I dx = M, 


5u dx = 0, 


(52) 


where is the support of u. Then to leading order, the difference in energy between 
the equilibrium solution and the perturbed one is 

AE = E{u + 5u) — E{u) K, / ■ 5udx. (53) 

In order for the energy to be stationary, AE must vanish for all perturbations 5u. 
As 5u has zero mass, a necessary and sufficient condition is that 6E/du in (53) is a 
constant, which we call A. That is. 


^ = A 
Su 


(54) 


in the support of u. The quantity A has a physical interpretation: it is the energy 
per unit mass at each point in space The condition ([54| is necessary for u to 
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Fig. 3. Numerically calculated energy minimizers in two dimensions. Density u is indicated 
by color. Mini mize rs of the nonlocal model appear in (a) and (c), and minimizers of the Cahn- 
Hilliard model l[29\ l appear in (b) and (d).(a,b) Mass M = 70, periodic box side length L = 40, 
calculated with N = 160 x 160 grid points. The minimizers have a well-defined peak. (c,d) Like 
(a,b) but M = 600. 


be a minimizer. For sufficient conditions, one must examine the continuation of A 
outside of the support as well as the second variation of the energy. These ideas are 
discussed at length in mi- While we have presented numerical results thus far, for 
the remainder of this section, we analyze the minimizers of (281 and (29). 


4.2. Minimizers of the nonlocal equation in one dimension. For (281, the 


energy (461 in one dimension becomes 


f 1 1 

E{u) = / -u[Qa*u]-\—+ Wudx. 
Jn‘2 6 


(55) 


From (47) and (54), and for W = 0, minimizers satisfy 


SE 

Su 



+ Qa * rt — A. 


(56) 


To make analytical progress, we consider the small and large mass limits in turn, with 
some ideas following [53J . 

When M is small, the characteristic length scale of the solution is small. Taylor 
expand the kernel as 




Q, 


(57) 
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Fig. 4. Radial profiles of th e ra dially-symmetric two dimensional minimizers i n F igure [3| 
Minimizers of the nonlocal model {ggp are dashed, and for the Cahn-Hilliard model solid, 

(a) Radial profiles for Figure\^a,b). (b) Radial profiles for Figure^^c,d). The plateau-like structure 
has peak density ||ii||oo ~ 1.5. 


Under this approximation, (56) becomes 


1 

2 


M 

Y 


1 

2 


\x\ * u = A, 


(58) 


where we have used the fact that l*u = M. Since (|a:|)a;a; = 2S{x), differentiate twice 
to obtain 


+ u = 0. (59) 

Without loss of generality, place the maximum of u at the origin; this implies Ux(0) = 0 

1 /2 

by symmetry. For convenience, rescale u by u(0) = ||u||oo and x by |jM||oo ■ Let 
P = u/\\u\\oo and ^ = a:/||u||oi , yielding 

(p')« + 2p = 0, p(0) = l, pc(0)=0. (60) 

Multiply by (p^)^, integrate once, and apply the conditions at ^ = 0 to obtain the 
constant of integration, leading to 





2 

3 



(61) 


Take the square root of both sides and separate variables to find an implicit solution. 





pdp 

y/l 


(62) 


As p —>■ 0, ^ approaches the half-width in our rescaled variable, which implies 



pdp 

\/l -P^ 



(63) 


II supp(u)|| = 2 ||m ||^2 
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The mass constraint is 

/'ll supp(«)||/2 

M = 2 


u{x) dx = (64) 


By combining (63l and (64), write || supp(m)|| and ||u||oo as functions of mass M, 

1/3 


fu oo = 


~ 0.721M2/3, 


supp(u)|| = \/-r - r - 


0 


1/6 


« 1.794Mi/^ 


(65a) 

(65b) 


For large mass, recall from Figuresandthat solutions approach wide plateaus 
with peak density ||m||oo ~ 1-5. To understand this behavior, we analyze the large M 
limit. When M is large, spatial scales are large. Thus Qa ~ ~<5(x) and (56) becomes 

5E 1 




( 66 ) 


from which it immediately follows that u is constant on the support of the solution. To 
determine what value that constant has, explicitly minimize the energy (46) subject 
to the mass constraint. The energy is 



supp(u)||. 


(67) 


However, the mass constraint dictates that for a this solution, ||supp(it)|| = M/u, 
which we substitute to obtain 


Thus 



( 68 ) 


(69) 


which yields the critical point u = 3/2. Since the solution is constant, ||m|||oo = 3/2, 
and II supp(u)|| = 2M/3. A curious feature of the approximation we have used is that 
it is invariant under any area-preserving map of the spatial domain; put differently, 
the approximation suggests that any arrangement of disjoint clumps having constant 
density u = 3/2 has the same energy. However, we know there is an energetic cost for 
each clump (neglected in our analysis) proportional to the measure of its perimeter. 
Thus, we expect that the global minimizer is a rectangular profile of density u = 3/2. 


4.3. Minimizers of the Cahn-Hilliard equation in one dimension. For 

the Cahn-Hilliard model (29), the energy (46) in one dimension becomes 


E{u) 



-I- -I- Wudx. 

2 6 


(70) 


As discussed in Section ( |^ has compactly supported solutions with zero contact 
slope. Here, we establish that energy minimizers share this property. We derive this 
result in one dimension, though one may extend the calculation to higher dimensions. 
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Consider an energy-minimizing solution u(x) with a support x G [a,/?]- We 
compute the first variation with respect to changes in the endpoints in addition to 
the solution itself. Let a = a + Sa,/3 = /3 + 5(3, u = u + 5u. Computing the first 
variation of the energy, we find 


6E = 


6E 
^ Su 


Su dx ■ 




6a, 


(71) 


where 


6E 1 2 
OU 2 


w. 


(72) 


In deriving we have used the fact that u vanishes at the endpoints x = a and 
X = (3. For SE to vanish, it is necessary that Ux{a) = Ux{i3) = 0. To see this, choose 
6u = 0 and suppose Ux(a) and Ux{P) are nonzero. In this case, choosing 5a > 0 would 
reduce the energy, as would choosing 5(3 < 0. (The change in mass resulting from this 
shifting of the endpoints slightly into the support is a second-order effect.) However, 
this violates the assumption that it is a minimizer. Therefore, Ux{a) = Ux{(3) = 0. 


Throughout much of the rest of our analysis of (29), we make use of the result we 


have just shown, namely that energy-minimizing solutions have zero contact slope. 


From (47) and (54), and for W = 0, minimizers satisfy 


5E 1 2 

— u — Uxx 

5u 2 


= A. 


(73) 


This differential equation has arisen previously in studies of the Cahn-Hilliard equa¬ 
tion [44] and the thin film equation [50] and has solutions that can be expressed in 
terms of elliptic functions. Multiply through by Ux, integrate once, and enforce that 
the energy minimizing solution has zero contact slope to obtain 


ul 

- -—-Au-^ = 0. 

6 2 2 


(74) 


By symmetry, UxiO) = 0 from which it follows that 


(75) 


at the center of the clump solution. Therefore, we know the relationship between 
||m||oo and A, namely 


A=>lL-2l 


(76) 


We derive expressions for the size of the support of the solution || supp(m)|| and for 
the mass M in terms of ||m||oo- For convenience, rearrange (74) as 


du\'^ 2 «^ 

- -2A„. 


(77) 
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For II supp(m)||, we have 


|isupp(M)|| 



supp(u)||/2 

dx, 



du 

dujdx' 



12 

V9-3||u|| 


K 


OO 


( I Ikiloo 

i,V3-ii«ii 


(78) 

(79) 

(80) 

(81) 


where K represents the complete elliptic integral of the first kind. For mass M we 
have 


I supp(ii)||/2 


M = 2 

= 2 

= 2 


udx, 

u du 
du/dx ’ 

udu 


f - _ 2Xu 


= V9-3||u|| 


I Halloo 
3- IImIIc 


-E 


( I Ikiloo \ 
3- Halloo y 


(82) 

(83) 

(84) 

(85) 


where E represents the complete elliptic integral of the second kind. As ||u||oo increases 
from 0 to 3/2, the mass M increases from 0 to infinity, as seen in Figure [^a). 

Ideally, we would invert these relationships to solve for ||u||oo and || supp(m)|| as a 
function of M, but this is cumbersome due to the elliptic integral functions. Therefore, 
as we did previously for the nonlocal model (28), we will analyze the small and large 
mass limits, which yield simple explicit expressions for ||u||oo and || supp(u)||. 

For small masses, neglect the nonlinear term in (73) to obtain 


U -\- Uxx — 


( 86 ) 


As mentioned above, the solution that minimizes E has zero contact angle and must 
also respect the mass constraint. Up to translation, the unique solution is 


u{x) 


^(1 + COSX) l^l < TT 

0 |a;| > 0 ■ 


(87) 


Thus, ||u||oo = M/n and || supp(u)|| = 27r. 

For large masses, the Uxx term in (73) is negligible and the solution is identical 
to the large mass limit of the nonlocal case considered above. That is, as M ^ oo. 
Halloo —>■ 3/2 and || supp(m)|| —)■ 2M/3. 

4.4. Minimizers in two dimensions. We briefly discuss the small and large 
mass limits in two dimensions. The small mass limit for the nonlocal model (28) 
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is difficult to compute explicitly because of convolution with the kernel Qa oc e 
An alternative approach is to use a different kernel such as the Bessel kernel used in 
PUSH- For our kernel, numerical explorations and a scaling argument both suggest 
that in the limit M —>■ 0, ||supp(u)|| and ||u||oo both approach zero in a fashion 
analogous to the problem in one dimension. The small mass limit for the Cahn-Hilliard 


model (29) is also analogous to its corresponding one dimensional case, satisfying 

u + = —A. 


The solution is 


L{r) = 


( M 

\ Mr) ' 
Mr*). 

< 

lo 



r <r^ 
r >r» 


( 88 ) 


(89) 


where the radius of the support, r*, is the first zero of Ji{r), and Jq and Ji are Bessel 
functions of the first kind. The support and peak density are 


I supp(m)|| = 7rr^, ||u||oo = 


M 


1 - 


1 


(90) 


Jo(r,)_ 

The large mass calculations for two dimensions are nearly identical to those for 

—>■3/2 and the minimizer is a 


one dimension. For 128^ and (291 as M —> oo 


disc of density u = 3/2, having radius 


r = 


2M 

37r 


(91) 


We highlight two key results from this section. First, the Cahn-Hilliard model 


(291 approximates the nonlocal model (281 well, especially for large masses and away 


from an aggregation’s boundaries. Second, in the limit of large mass, both models 
produce groups with an energetically preferred peak density of 3/2. This value is 
the benchmark against which we address a biologically motivated question: can an 
external potential modeling the environment be used to reduce the peak density of an 
aggregation? 


5. Energy minimizers in an external potential well. With (291 established 


as a reasonable approximation of (28), at least for large enough masses and away 


from transition regions, we now use the degenerate Cahn-Hilliard model as a testbed 
to investigate the effects of the environment; that is, we now allow W ^ 0. More 
specifically, we focus on cases where VF is a potential well, and show that when this 
well has sufficiently tall and steep sides, it can be replaced by an approximate obstacle 
potential that is either zero or a sufficiently large, and in fact effectively infinite, 
constant. In this scenario, the support of a minimizer is constrained to avoid locations 
where W is large. For the remainder of this paper, we simply refer to “obstacle 
potentials,” but remind the reader that in our computational implementation, these 
are approximate obstacle potentials with large, rather than infinite, height. 


The velocity in (29) only depends on the gradient of W. Furthermore, adding a 
constant to W affects the energy (70) by adding a multiple of the conserved mass. 
Therefore, W is arbitrary up to a constant. Without loss of generality, take the 
minimum value of W to be zero. As a prototypical example, consider a Gaussian 
barrier placed at each end of the interval [0,1/], that is, 

A 


W{x) = 


2i/7ro 




-{L-xYHa^ 


(92) 
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where A measures the volume of the barriers and a is the characteristic width. 

Figure shows the ensuing minimizers of (291. Panel (a) varies the height A 
with fixed width a = L/40 on a domain of length L = 10 with mass M = 15. As A 
increases, the density at the boundary decreases and is eventually driven to rt = 0. 
For sufficiently large A, W drives the minimizer to zero over an interval, creating 
vacancies near the barriers. Panel (b) is similar, but varies a with fixed A = 10, 
which is sufficiently large to ensure u = 0 at the boundary. As the barrier width a 
decreases, the vacant interval shrinks to a point, and W pins the minimizer to u = 0 
at the boundaries. For the smallest few values of a, W is nonzero only at a the edges 
of the computational domain; that is to say, W approaches infinity at the domain 
endpoints and is zero within numerical tolerance elsewhere. 

To restate: as we decrease cr, the narrow Gaussian, which we have numerically 
under-resolved, approaches an obstacle potential, and the minimizer approaches one 
which is pinned to m = 0 at the boundary. This limit suggests that we may replace 
our smooth potential W with an approximate obstacle potential that is large on some 
subset of the gridpoints and zero elsewhere. In this limit, continuity of W is lost. One 
consequence is that the analysis in Section |4.3| demonstrating zero contact slope for 
minimizers no longer holds. Indeed, it is apparent in Figurej^b) that minimizers may 
have nonzero slope where they are pinned by the potential at the domain boundary. 

The strategy of adopting an obstacle potential is effective in two dimensions as 
well. Figure]^ shows a minimizer of (29) with mass M = 9 on a domain with sides of 
length L = 3. In the absence of an external potential, the minimizer has constant den¬ 
sity u(x) = 1. Panel (a) takes the potential kF(x) = 10^ at computational gridpoints 
that trace out the white heart-shaped curve, and W = 0 elsewhere. The potential 
drives the minimizer to m = 0 along this curve, apparent as the dark blue region, 
and the peak density is ||w||oo ~ 2.2. Panel (b) is similar, but takes fF(x) = 10^ 
on and exterior to the curve, creating a heart-shaped well. The potential confines 
the minimizer to the well, drives it to zero along the boundary, and pushes the peak 
density even higher, to ||u||oo ~ 4.7. 

The next section focuses on cases where W is an obstacle potential, and examines 
strategies for reducing the peak density by engineering W. 


6 . The effect on peak density of random versus periodic potentials. Our 

biological motivation for this study is the dynamics of locust populations. Intuitively, 
one might think that the distribution of many organisms in the wild is affected by 
the distribution of resources on which they subsist, and indeed, this issue is part 
of the rich field of biogeography. The way that resources shape the distribution of 
locusts is intimately tied to locusts’ phase polyphenism. Phase polyphenism refers 
to individuals of the same genotype manifesting different phenotypes [U I60j . For 
locusts, a pivotal phenotypic difference is the insect’s social behavior, or lack thereof. 
A locust may exist in a solitarious phase in which it seeks isolation or a gregarious 
phase in which it seeks other locusts. Locusts in the gregarious phase may form large, 
migratory swarms that decimate crops. Therefore, dense social aggregations on the 
ground are of interest as they are precursors to dangerous flying groups. 

Previous work identifies a population density threshold above which a collec¬ 
tive transition to a dangerous, gregarious group takes place. Thus, one way to prevent 
airborne swarms might be to suppress the density of locusts on the ground, so that 
the group never reaches the critical threshold. Because social locusts display taxis not 
only to each other, but to food, one can ask whether there exist external potentials 
W, modeling food sources, which serve to minimize the peak density ||m||oo of the 
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Fig. 5. Minimizers of the degenerate Cahn-Hilliard equation with Gaussian potential 
barriers at the domain edges, as described by W{x) in l[92f . Domain length L = 10 with N = 160 
computational grid points. Total mass is M = 15. (a) Fixed Gaussian width a = L/40 for potential 
volume A varying between 1 and 10. As the amplitude increases, ii(0) and u(L) are driven to zero 
and for sufficiently large A there is a vacant segment near the potential barrier, (b) Fixed Gaussian 
amplitude A = 10 for width a varying logarithmically between 0.0025 and 0.25. Here, A is sufficiently 
large to drive u to zero near the boundary. As a decreases, the vacant segment shrinks. The smallest 
few values of a to yield effectively nonzero W only at a the edges of the computational domain; that 
is to say, W is large at the domain endpoints and zero within numerical tolerance elsewhere. This 
choice pins the minimizer to be zero just at the endpoints, where a finite contact slope is observed. 


population. We investigate this problem within the framework of (29). The locust 
model in [5^ includes additional effects - crucially, the phase change from solitarious 
to gregarious locusts - but an understanding of the interaction between locusts and 
food sources even in the absence of phase change is an appropriate starting point for 
investigation. 

Figure reveals some of the intricacies of the problem. Consider a one dimen¬ 
sional domain and obstacle potentials W (x) describing a sequence of square wells with 
dividers between them, intended to model cropland divided into plots. Each panel 
shows a minimizer of (291. At low and moderate masses M (top and middle rows), 
applying a large number of square wells in W (x) to the fixed domain reduces the 
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Fig. 6. Minimizers of the degenerate Cahn-Hilliard equation for two obstacle potentials 
VF(x). The domain has sides of length L = 3 with N = 160 computational grid points along each 
axis. Total mass is M = 9. (a) lF(x) = 10^ on the heart-shaped curve, and W = 0 elsewhere. The 
potential pins the minimizers to u = 0 along this curve, (h) Similar to (a), but the potential is now 
nonzero on and exterior to the curve, creating a heart-shaped well. In the absence of a potential, 
the minimizer is a constant, u{x) = 1. In (a), W drives the peak density higher, to ||w||cx 3 ~ 2.2. In 
(b), W drives the peak density even higher, to Halloo ~ 4.7. 


peak density ||u||oo, whereas at high masses (bottom row) it has the opposite effect. 
Furthermore, at low and moderate densities, social aggregation can cause the mass 
to clump in a subset of the wells, as seen in panels (b) and (e), which may serve to 
either reduce or increase peak density. 

In summary, the goal of this section is to choose external potentials W to min¬ 
imize the peak density ||u||oo using the Cahn-Hilliard model (29) as a testbed. For 


tractability, we restrict attention to the class of obstacle potentials. The next three 
subsections focus on a one dimensional domain with periodic W. In the the final sub¬ 
section, we compare results for periodic potentials with those for random potentials 
in both one and two dimensions. Random potentials tend to increase peak density 
||m||oo while periodic potentials can decrease it. 

6.1. Energy minimizers in a single square well. For the Cahn-Hilliard 
model (29), when ||supp(u)|| for a given mass on an infinite domain is less than 
the well width £, the minimizer is identical to that found in Section However, 
when the support is larger, the minimizer feels the boundaries of the square well. For 
sufficiently large W at the edge of the well (that is, sufficient contrast between desert 
and planted regions), u will be driven to zero. 


From (47) and (|54[) , minimizers satisfy 


5E 


1 


— 'tt U 'ttxx 

5u 2 


= A, 


(93) 


subject to the mass constraint and the added restriction that u = 0 at the edges of 
the well because of the external potential. An implicit solution for m as a function of 
A, mass M, and well size ^ could be found in terms of elliptic integrals. However it 
is simpler to consider small and large mass limits, as in Section]^ and augment our 
calculations with numerical computations. First, we calculate the energy E of the 
minimizing solution for several cases depending on M and 1. 
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Fig. 7. Minimizers of the Cahn-Hilliard model f29]) on a domain of length L = 127r for varying 
mass M and varying numbers of square wells in the external potential W{x). (a-c) M = ir. (d-f) 
M = An. (g-i) M = IGtt. In the first column, there is one square well, an arrangement we use in 
analogy to an isolated plot of planted land. In the second column, there are eight square wells, and 
in the third column there are 24- At low and intermediate M (top and middle rows), applying a 
large number of square wells in W{x) to the fixed domain reduces the peak density ||ii||ooj o-s seen 
in (c) and (f). In contrast, at large M (bottom row), increasing the number of square wells applied 
via W{x) increases ||w||oo- Eor the numerical computations we use N = 744 grid points. 


Case lA: Small M, £ < 27r. In Section]^ we showed that for an infinite domain, 
the small M solutions have support of size 27r. We expect that for our present case, 
the solution will occupy the entire domain. To examine the small mass limit, linearize 
(931 to obtain the boundary value problem 

U Uxx — 


u(0) = u{£) = 0, 

in addition to the mass constraint. The solution is 

^cos(a: — ^/2) — cos(£/2) 


(94) 


u{x) = M- 


2 sin(£/2) — £ cos(^/2) 


(95) 


Direct computation via substitution into (701 yields the energy for this solution, 

M2 


E = 


cos(f/2) 

~2~ 2sin(G2 )-£cos(£/2) 


+ C>(M^). 


(96) 


Crucially, the energy changes from positive to negative as £ increases through tt. 
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Case IB: Small M, i > 27r. The results are identical to the infinite interval result 
(so long as the support of the solution is less than t). Substituting (87) into (70), we 
find the energy 


E = - 


47r 




(97) 


Case II: Large M. Numerical simulation shows the density profile to be constant 
in the potential well except for a small boundary layer near the sides of the well of 
thickness 0{\fljM). Exploiting this fact, we have 

u{x) = Mji, (98) 


and thus the energy is 


E = 


m3 

W 


O 


/M5/2\ 

) ■ 


(99) 


6.2. Energy per unit mass in a single well. Momentarily, we will make 
predictions about energy minimizers in multiple wells. To minimize the total energy 
we will show that a system with multiple wells distributes mass among those wells to 
minimize the energy per unit mass. Thus, as another preliminary step, we use the 
results of the previous subsection to calculate the energy per unit mass in a single 
well. Define the energy per unit mass, 


£{M) = E{M)/M. 


( 100 ) 


From the previous subsection, we can determine the behavior of £{M) for the limiting 
cases of small and large M. We will see that £{M) may have a minimum, which will 
be a crucial characteristic when we consider the case of multiple wells. We refer to 
the value of M that minimizes £{M) as M*. 

Case I: £ < TT. A sample computation of £{M) is shown in Figure [^a). It is 
monotonically increasing. This result is consistent with the small mass result ( |96[ ), 
from which it follows that 


£(M) = 


M 


cos(£/2) 


2 2sin(£/2) — £cos(£/2) ’ 


which is linearly increasing with M. Also, from the large mass result (99), 


£{M) = 


M2 

6^’ 


( 101 ) 


( 102 ) 


which is again increasing. We conclude that M* = 0, that is, having the minimum 
possible mass in the well minimizes £{M). 

Case II: TT < £ < 27r. A sample calculation of £{M) is shown in Figure |^d). 
Relations (1011 and (102) still hold, but crucially, the coefficient in the small mass 
limit is negative rather than positive, as it was in Case I. Thus, we expect M* to 
occur between the decreasing behavior at small mass and the increasing behavior at 
large mass, as evidenced in the figure. 

Case III: £ > 27r (but finite). Recall that for an infinite domain and in the limit of 


small mass, the minimizer is given by (87) and the support is 27r. It is not surprising. 
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Fig. 8. Properties of minimizers of the Cahn-Hilliard model ^29^ with a single square well of 
width i used for the external potential W{x). Columns are the energy per unit mass S(M), the peak 
density ||ti||oo? o-nd the size of the support || supp(ti)||, all as a function of mass M. (a-c) Case I, 
i = 7r/2. Here, S{M) and ||ii||oo increases monotonically and supp(u) fills the entire domain, (d-f) 
Case II, i = 37r/2. Here, S{M) reaches a minimum at M = M*, ||'ii||oo increases monotonically 
and supp(ii) fills the entire domain, (g-i) Case III, i = 57r/2. Here, S{M) reaches a minimum at 
M = M*, ||w||oo increase monotonically and supp{u) widens and eventually fills the entire domain, 
(j-l) Case IV, i chosen to be effectively infinite, much larger than supp(u). Here, S{M) decreases 
monotonically to —3/8 and ||w||oo increases monotonically to 3/2. Additionally, supp(ti) widens 
without bound and, asymptotically, grows linearly. 


therefore, that the same small mass result holds for intervals wider than 27r. Plugging 
the minimizer into the energy and dividing by M yields 


£{M) = - 


M 

47r 


(103) 
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which is decreasing. As the mass increases, the support widens as well, and eventually 
fills the domain, as shown in Figure |^i). For sufficiently large mass, the density is 
nearly constant with u{x) « M/£ (except at the boundaries) and the approximation 


(1021 still holds. As in Case II above, E{M) is decreasing at small mass, increasing 
at large mass and reaches a minimum at M* between these two regimes. 


Case IV: f large (effectively infinite). The small mass approximation (1031 still 


holds. As M increases, the minimizer follows the branch of solutions found for the infi¬ 
nite interval. Moreover, by the energy argument in Sectionj^ the solution approaches 
a single rectangle of height u = 3/2, having, by direct calculation. 


S{M) = 


(104) 


Thus, £{M) decreases monotonically but approaches this constant, and so M* is 
effectively infinite. It will be energetically favorable for the mass to form a clump of 
density u = 3/2 inside a single well as was seen previously. 

6.3. Energy minimizers in multiple -wells. We now take the external poten¬ 
tial W to consist of a periodic sequence of square wells of width £. The key question 
is whether the mass will distribute itself equally amongst the wells or whether it is 
energetically preferred to concentrate in a subset of the wells. We use the results from 
the previous subsections to answer this question. 

More concretely, suppose that there are n equal wells in the domain and that for 
each well, E{m) is the energy associated with a minimizer of mass m located in that 
well. The total energy for a distribution of masses rui > 0 is 

n 

E{M) = Y,E{mi), (105) 

where 

n 

(106) 

i=l 

because total mass is fixed. To proceed, consider an ansatz where the mass is equally 
distributed between k of the n wells. In this case, 

E{M) = kE[M/k) = —E{fi) = (107) 

where fj, is the mass contained in each occupied well. From this equation, we deduce 
that a necessary condition to minimize E{M) is to minimize the energy per unit mass 
£’(/r) over admitted values of p,. 

For the sake of argument, suppose one can choose ^ = M*. In this case, E{M) is 
globally minimized. In practice, the values of /r are quantized since fi = M/k where 
fc is a positive integer less than or equal to n. For large enough n, the system appears 
to choose /X that approximates M* well. 

Now apply these ideas to the numerical results in Figure For panels (a,d,g) in 
the first column, there is one square well, that is, £ = L = 127r. For panels (b,e,h) in 
the second column, there are eight square wells, that is £ = L/8 = 37r/2, for which we 
have estimated M* « 2.97 from the data in Figure |^d). In panel (b), as M « M*, 
the mass fills a single well. In panel (e), as AM* < M < 5M*, the mass fills five 
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of the wells. In panel (h), as M > 8M*, the mass equipartitions among all eight 
wells. Finally, for panels (c,f,i) in the third column, there are 24 square wells, that is, 
£ = L/24 = 7r/2 for which M* = 0. Therefore, for all values of M, mass equipartitions 
among all wells. 

6.4. Comparing random and periodic potentials. We now demonstrate 
that a periodic pattern of resources in the environment can reduce peak density below 
what might be seen with a random distribution of resources. This effect is evident 
in Figure]^ which compares these two cases in one and two dimensions. For all four 
panels, the average density is u = 1/2, which is below the peak density of 3/2 seen in 
minimizers in the absence of an external potential (see Section]^. Additionally, we fix 
(j), the volume fraction of the computational domain on which the obstacle potential 
W is effectively infinite (here, 10^) and we take W = 0 elsewhere. These choices of 
W model, respectively, regions barren and lush with vegetation. In Figure]^ we set 

0.08. 

First, consider the one dimensional case. In panel (a), mass equipartitions among 
all 40 available potential wells, and ||m||oo ~ 0.76. For the randomly distributed 
potential used in (b), some wells are left vacant and the peak density increases to 
||w||oo ~ 1.95. Now, for two dimensions, panel (c) demonstrates that mass equiparti¬ 
tions among all 25 potential wells in the five-by-five grid, and ||m||oo ~ 1-2. However, 
for the randomly distributed potential used in (d), there exist hot spots with peak 
densities as high as ||u||oo ~ 1-7- 

In short, periodic potentials can suppress the peak density below 3/2, while ran¬ 
dom potentials can concentrate density at values greater than 3/2. Figure [l0| provides 
a systematic investigation for varying volume fraction in one and two dimensions. 
In both cases, at low cj), ||m||oo ~ 3/2 (dotted line) for both random (shown as squares) 
and periodic (shown as circles) potentials. As (j) increases, the periodic potential sup¬ 
presses peak density below 3/2 whereas the random potential increases it above. This 
result suggests the viability of carefully engineering resource layout as a strategy for 
reducing a population’s peak density. 

7. Conclusion. Understanding how populations respond to social forces and 
environmental cues such as resource distribution is essential for modeling biological 
groups. For example, for desert locusts, spatially varying resource distributions can 
concentrate populations, which in the wild can lead to collective gregarization and 
dangerous outbreaks [36] . We have laid out a framework for studying these problems 
by exploiting energy minimization of the governing models. 

More specifically, we began with aggregation model describing nonlocal social 
attraction, local, nonlinear repulsion, and an external potential modeling the environ¬ 
ment. This model possesses an energy that is minimized by the dynamics. However, 
the model poses (at least) three challenges. First, time-dependent simulations are 
costly, and may converge slowly to attractors. Second, attractors typically have com¬ 
pact support. Standard simulation methods may produce oscillations emanating from 
contact points which typically violates nonnegativity of the solution. Third, the non¬ 
locality couples all spatial locations, rendering numerical computations expensive, 
especially in two or more dimensions. 

To address the first two challenges, we exploit the energy minimizing dynamics 
by focusing on the analysis and numerical computation of minimizers. This avoids 
costly computations of dynamics, and additionally, minimization algorithms allow us 
to enforce nonnegativity of the solution as a simple constraint in the optimization 
procedure. To address the third challenge, we perform a long-wave approximation of 
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Fig. 9. Example minimizers of ^29^ for random and periodic obstacle potentials iy(x). For 
these potentials, we fix 4>, the volume fraction of the computational domain in which the environ¬ 
mental potential W is effectively infinite (here, Elsewhere, W = 0. These choices of W model, 

respectively, regions barren and lush with vegetation. (a,b) Minimizers in one dimension for M = 50 
on a domain of length L = 100 with 480 computational gridpoints and 0 = 1/12 ps 0.083. For the 
periodic potential used in (a), mass equipartitions into all fO potential wells and ||w||oo ~ 0.76. For 
the randomly distributed potential used in (b), some wells are left vacant and the peak density in¬ 
creases to ||w||oo ~ 1.95. (c,d) Minimizers in two dimensions for M = 648 on a domain with sides of 
length L = 36 with N = 120 computational grid points along each axis and 0 = A1 Pi 0.082. For 
the periodic potential used in (c), mass equipartitions into all 25 potential wells in the five-by-five 
grid, and Halloo ~ 1-2. For the randomly distributed potential used in (d), there exist hot spots with 
peak densities as high as Halloo ~ 1-7. 


our model to obtain a fourth-order degenerate Cahn-Hilliard equation which is local 
in space while preserving the energy minimizing character of the original model. In 
summary, we have developed an approach to studying aggregations that replaces sim¬ 
ulating a nonlocal equation with minimizing a local energy subject to a nonnegativity 
constraint. 

In the absence of an external (that is, environmental) potential, minimizers of the 
nonlocal and Cahn-Hilliard equations agree well at large masses and away from the 
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<t> 



<t> 


Fig. 10. Peak density ||'u||oo of minimizers of the Cahn-Hilliard equation for two classes 
of obstacle potentials VF(x). We vary (j), the volume fraction of the computational domain in which 
the environmental potential W is effectively infinite (here, 10“^/ Elsewhere, W = 0. These choices of 
W model, respectively, regions barren and lush with vegetation. Open circles correspond to choosing 
for W a periodic pattern of square wells, representing engineered cropland with barriers between crop 
beds. Squares correspond to choosing for W a uniformly random distributed pattern, producing an 
irregular, clumpy vegetation landscape. Error bars represent maximum and minimum values over 
ten realizations of W for the random case. Dotted lines represent HiiHoo = 3/2, the energetically 
preferred peak density for large mass solutions in a large domain, (a) Minimizers in one dimension, 
with domain length L = 100, population mass M = 50, simulated with N = 480 gridpoints. (b) 
Minimizers in two dimensions, with domain length L = 36, population mass M = 648, and N = 120 
gridpoints along each axis. For both (a) and (b), the average population density is 1/2. In both 
panels, a randomly generated W tends to increase peak density ||t4||oo above 3/2 while periodic 
square wells tend to decrease it below 3/2. 


contact points. In the limit of large domains and modest densities, these minimizers 
approach a compactly supported clump of population with peak density 3/2 and steep 
edges. 

An external potential drives mass to accumulate near potential minima, often 
creating collections of compactly supported clumps. Many environmental landscapes 
can be partitioned into regions that are either inviting or inhospitable to biological 
organisms. This characterization arises naturally in numerical computations for suffi¬ 
ciently steep potential wells. Thus, we model these landscapes via obstacle potentials 
that are effectively infinite (inhospitable) on some subset of the domain and zero 
(inviting) elsewhere. 

In this obstacle potential testbed, we ask how resource distribution effects peak 
population density. In the absence of an external potential and at modest average 
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population densities, aggregations form with peak density 3/2, as mentioned above. 
A spatially random obstacle potential modeling the natural environment drives peak 
densities higher, even when the inhospitable portion of the domain is less than 10% 
of the total. However, a periodic pattern of square wells modeling divided crop beds 
can reduce the peak density below 3/2. 

Anthropologists have chronicled the use of periodic planting patterns since an¬ 
cient times, and typically attribute this strategy to better water management and 
thus increased crop yield. We are intrigued, but admit that it is pure speculation, 
to ask whether this strategy might confer any pest control advantages. Regardless, 
the biological literature notes that for locusts, “fractal dimension” of resources in the 
environment can drive large peak population densities [36]. While perhaps just a 
caricature, our simple model suggests that carefully distributing resources can inhibit 
populations from reaching their intrinsic peak density and plausibly could help con¬ 
trol undesirable biological phenomena such as gregarious locust outbreaks that are 
triggered by surpassing a density threshold. 
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